void shake3(celldata_t *cell, int i, real dtfsq, real tol, int maxiter){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  int ig2 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[1]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  real rmass2 = cell->rmass[ig2];
  real d01sq = cell->shake[i].rsq[0];
  real d02sq = cell->shake[i].rsq[1];
  vec<real> r01, ruc01, r02, ruc02;
  r01 = cell->x[ig0] - cell->x[ig1];
  r02 = cell->x[ig0] - cell->x[ig2];
  ruc01 = cell->shake_xuc[ig0] - cell->shake_xuc[ig1];
  ruc02 = cell->shake_xuc[ig0] - cell->shake_xuc[ig2];
  real ruc01sq = ruc01.norm();
  real ruc02sq = ruc02.norm();
  real r0101 = r01.dot(r01);
  real r0102 = r01.dot(r02);
  real r0202 = r02.dot(r02);
  real a00 = 2 * r01.dot(ruc01)*(+rmass0+rmass1);
  real a01 = 2 * r01.dot(ruc02)*(+rmass0);
  real a10 = 2 * r02.dot(ruc01)*(+rmass0);
  real a11 = 2 * r02.dot(ruc02)*(+rmass0+rmass2);
  real invmatdet = 1/(a00*a11 - a01*a10);
  real a00inv = (a11) * invmatdet;
  real a01inv = (-a01) * invmatdet;
  real a10inv = (-a10) * invmatdet;
  real a11inv = (a00) * invmatdet;
  real quad0_00 = r0101*(rmass0 + rmass1)*(rmass0 + rmass1);
  real quad0_01 = 2*r0102*rmass0*(rmass0 + rmass1);
  real quad0_11 = r0202*(rmass0)*(rmass0);
  real quad1_00 = r0101*(rmass0)*(rmass0);
  real quad1_01 = 2*r0102*rmass0*(rmass0 + rmass2);
  real quad1_11 = r0202*(rmass0 + rmass2)*(rmass0 + rmass2);
  real l0 = 0;
  real l1 = 0;
  int done = 0;
  for (int i = 0; i < maxiter; i ++) {
    real quad0 = quad0_00*l0*l0 + quad0_01*l0*l1 + quad0_11*l1*l1;
    real quad1 = quad1_00*l0*l0 + quad1_01*l0*l1 + quad1_11*l1*l1;
    real b0 = d01sq - ruc01sq - quad0;
    real b1 = d02sq - ruc02sq - quad1;
    real l0_new = a00inv*b0 + a01inv*b1;
    real l1_new = a10inv*b0 + a11inv*b1;
    int almost_done = fabs(l0_new - l0) < tol && fabs(l1_new - l1) < tol;
    l0 = l0_new;
    l1 = l1_new;
    if (done) break;
    done = almost_done;
  }
  l0 = l0 / dtfsq;
  l1 = l1 / dtfsq;
  cell->f[ig0] += r01 * l0;
  cell->f[ig0] += r02 * l1;
  cell->f[ig1] -= r01 * l0;
  cell->f[ig2] -= r02 * l1;
}
void shake3angle(celldata_t *cell, int i, real dtfsq, real tol, int maxiter){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  int ig2 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[1]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  real rmass2 = cell->rmass[ig2];
  real d01sq = cell->shake[i].rsq[0];
  real d02sq = cell->shake[i].rsq[1];
  real d12sq = cell->shake[i].rsq[2];
  vec<real> r01, ruc01, r02, ruc02, r12, ruc12;
  r01 = cell->x[ig0] - cell->x[ig1];
  r02 = cell->x[ig0] - cell->x[ig2];
  r12 = cell->x[ig1] - cell->x[ig2];
  ruc01 = cell->shake_xuc[ig0] - cell->shake_xuc[ig1];
  ruc02 = cell->shake_xuc[ig0] - cell->shake_xuc[ig2];
  ruc12 = cell->shake_xuc[ig1] - cell->shake_xuc[ig2];
  real ruc01sq = ruc01.norm();
  real ruc02sq = ruc02.norm();
  real ruc12sq = ruc12.norm();
  real r0101 = r01.dot(r01);
  real r0102 = r01.dot(r02);
  real r0112 = r01.dot(r12);
  real r0202 = r02.dot(r02);
  real r0212 = r02.dot(r12);
  real r1212 = r12.dot(r12);
  real a00 = 2 * r01.dot(ruc01)*(+rmass0+rmass1);
  real a01 = 2 * r01.dot(ruc02)*(+rmass0);
  real a02 = 2 * r01.dot(ruc12)*(-rmass1);
  real a10 = 2 * r02.dot(ruc01)*(+rmass0);
  real a11 = 2 * r02.dot(ruc02)*(+rmass0+rmass2);
  real a12 = 2 * r02.dot(ruc12)*(+rmass2);
  real a20 = 2 * r12.dot(ruc01)*(-rmass1);
  real a21 = 2 * r12.dot(ruc02)*(+rmass2);
  real a22 = 2 * r12.dot(ruc12)*(+rmass1+rmass2);
  real invmatdet = 1/(a00*a11*a22 - a00*a12*a21 - a01*a10*a22 + a01*a12*a20 + a02*a10*a21 - a02*a11*a20);
  real a00inv = (a11*a22 - a12*a21) * invmatdet;
  real a01inv = (-a01*a22 + a02*a21) * invmatdet;
  real a02inv = (a01*a12 - a02*a11) * invmatdet;
  real a10inv = (-a10*a22 + a12*a20) * invmatdet;
  real a11inv = (a00*a22 - a02*a20) * invmatdet;
  real a12inv = (-a00*a12 + a02*a10) * invmatdet;
  real a20inv = (a10*a21 - a11*a20) * invmatdet;
  real a21inv = (-a00*a21 + a01*a20) * invmatdet;
  real a22inv = (a00*a11 - a01*a10) * invmatdet;
  real quad0_00 = r0101*(rmass0 + rmass1)*(rmass0 + rmass1);
  real quad0_01 = 2*r0102*rmass0*(rmass0 + rmass1);
  real quad0_02 = -2*r0112*rmass1*(rmass0 + rmass1);
  real quad0_11 = r0202*(rmass0)*(rmass0);
  real quad0_12 = -2*r0212*rmass0*rmass1;
  real quad0_22 = r1212*(rmass1)*(rmass1);
  real quad1_00 = r0101*(rmass0)*(rmass0);
  real quad1_01 = 2*r0102*rmass0*(rmass0 + rmass2);
  real quad1_02 = 2*r0112*rmass0*rmass2;
  real quad1_11 = r0202*(rmass0 + rmass2)*(rmass0 + rmass2);
  real quad1_12 = 2*r0212*rmass2*(rmass0 + rmass2);
  real quad1_22 = r1212*(rmass2)*(rmass2);
  real quad2_00 = r0101*(rmass1)*(rmass1);
  real quad2_01 = -2*r0102*rmass1*rmass2;
  real quad2_02 = -2*r0112*rmass1*(rmass1 + rmass2);
  real quad2_11 = r0202*(rmass2)*(rmass2);
  real quad2_12 = 2*r0212*rmass2*(rmass1 + rmass2);
  real quad2_22 = r1212*(rmass1 + rmass2)*(rmass1 + rmass2);
  real l0 = 0;
  real l1 = 0;
  real l2 = 0;
  int done = 0;
  for (int i = 0; i < maxiter; i ++) {
    real quad0 = quad0_00*l0*l0 + quad0_01*l0*l1 + quad0_02*l0*l2 + quad0_11*l1*l1 + quad0_12*l1*l2 + quad0_22*l2*l2;
    real quad1 = quad1_00*l0*l0 + quad1_01*l0*l1 + quad1_02*l0*l2 + quad1_11*l1*l1 + quad1_12*l1*l2 + quad1_22*l2*l2;
    real quad2 = quad2_00*l0*l0 + quad2_01*l0*l1 + quad2_02*l0*l2 + quad2_11*l1*l1 + quad2_12*l1*l2 + quad2_22*l2*l2;
    real b0 = d01sq - ruc01sq - quad0;
    real b1 = d02sq - ruc02sq - quad1;
    real b2 = d12sq - ruc12sq - quad2;
    real l0_new = a00inv*b0 + a01inv*b1 + a02inv*b2;
    real l1_new = a10inv*b0 + a11inv*b1 + a12inv*b2;
    real l2_new = a20inv*b0 + a21inv*b1 + a22inv*b2;
    int almost_done = fabs(l0_new - l0) < tol && fabs(l1_new - l1) < tol && fabs(l2_new - l2) < tol;
    l0 = l0_new;
    l1 = l1_new;
    l2 = l2_new;
    if (done) break;
    done = almost_done;
  }
  l0 = l0 / dtfsq;
  l1 = l1 / dtfsq;
  l2 = l2 / dtfsq;
  cell->f[ig0] += r01 * l0;
  cell->f[ig0] += r02 * l1;
  cell->f[ig1] -= r01 * l0;
  cell->f[ig1] += r12 * l2;
  cell->f[ig2] -= r02 * l1;
  cell->f[ig2] -= r12 * l2;
}
void shake4(celldata_t *cell, int i, real dtfsq, real tol, int maxiter){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  int ig2 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[1]];
  int ig3 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[2]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  real rmass2 = cell->rmass[ig2];
  real rmass3 = cell->rmass[ig3];
  real d01sq = cell->shake[i].rsq[0];
  real d02sq = cell->shake[i].rsq[1];
  real d03sq = cell->shake[i].rsq[2];
  vec<real> r01, ruc01, r02, ruc02, r03, ruc03;
  r01 = cell->x[ig0] - cell->x[ig1];
  r02 = cell->x[ig0] - cell->x[ig2];
  r03 = cell->x[ig0] - cell->x[ig3];
  ruc01 = cell->shake_xuc[ig0] - cell->shake_xuc[ig1];
  ruc02 = cell->shake_xuc[ig0] - cell->shake_xuc[ig2];
  ruc03 = cell->shake_xuc[ig0] - cell->shake_xuc[ig3];
  real ruc01sq = ruc01.norm();
  real ruc02sq = ruc02.norm();
  real ruc03sq = ruc03.norm();
  real r0101 = r01.dot(r01);
  real r0102 = r01.dot(r02);
  real r0103 = r01.dot(r03);
  real r0202 = r02.dot(r02);
  real r0203 = r02.dot(r03);
  real r0303 = r03.dot(r03);
  real a00 = 2 * r01.dot(ruc01)*(+rmass0+rmass1);
  real a01 = 2 * r01.dot(ruc02)*(+rmass0);
  real a02 = 2 * r01.dot(ruc03)*(+rmass0);
  real a10 = 2 * r02.dot(ruc01)*(+rmass0);
  real a11 = 2 * r02.dot(ruc02)*(+rmass0+rmass2);
  real a12 = 2 * r02.dot(ruc03)*(+rmass0);
  real a20 = 2 * r03.dot(ruc01)*(+rmass0);
  real a21 = 2 * r03.dot(ruc02)*(+rmass0);
  real a22 = 2 * r03.dot(ruc03)*(+rmass0+rmass3);
  real invmatdet = 1/(a00*a11*a22 - a00*a12*a21 - a01*a10*a22 + a01*a12*a20 + a02*a10*a21 - a02*a11*a20);
  real a00inv = (a11*a22 - a12*a21) * invmatdet;
  real a01inv = (-a01*a22 + a02*a21) * invmatdet;
  real a02inv = (a01*a12 - a02*a11) * invmatdet;
  real a10inv = (-a10*a22 + a12*a20) * invmatdet;
  real a11inv = (a00*a22 - a02*a20) * invmatdet;
  real a12inv = (-a00*a12 + a02*a10) * invmatdet;
  real a20inv = (a10*a21 - a11*a20) * invmatdet;
  real a21inv = (-a00*a21 + a01*a20) * invmatdet;
  real a22inv = (a00*a11 - a01*a10) * invmatdet;
  real quad0_00 = r0101*(rmass0 + rmass1)*(rmass0 + rmass1);
  real quad0_01 = 2*r0102*rmass0*(rmass0 + rmass1);
  real quad0_02 = 2*r0103*rmass0*(rmass0 + rmass1);
  real quad0_11 = r0202*(rmass0)*(rmass0);
  real quad0_12 = 2*r0203*(rmass0)*(rmass0);
  real quad0_22 = r0303*(rmass0)*(rmass0);
  real quad1_00 = r0101*(rmass0)*(rmass0);
  real quad1_01 = 2*r0102*rmass0*(rmass0 + rmass2);
  real quad1_02 = 2*r0103*(rmass0)*(rmass0);
  real quad1_11 = r0202*(rmass0 + rmass2)*(rmass0 + rmass2);
  real quad1_12 = 2*r0203*rmass0*(rmass0 + rmass2);
  real quad1_22 = r0303*(rmass0)*(rmass0);
  real quad2_00 = r0101*(rmass0)*(rmass0);
  real quad2_01 = 2*r0102*(rmass0)*(rmass0);
  real quad2_02 = 2*r0103*rmass0*(rmass0 + rmass3);
  real quad2_11 = r0202*(rmass0)*(rmass0);
  real quad2_12 = 2*r0203*rmass0*(rmass0 + rmass3);
  real quad2_22 = r0303*(rmass0 + rmass3)*(rmass0 + rmass3);
  real l0 = 0;
  real l1 = 0;
  real l2 = 0;
  int done = 0;
  for (int i = 0; i < maxiter; i ++) {
    real quad0 = quad0_00*l0*l0 + quad0_01*l0*l1 + quad0_02*l0*l2 + quad0_11*l1*l1 + quad0_12*l1*l2 + quad0_22*l2*l2;
    real quad1 = quad1_00*l0*l0 + quad1_01*l0*l1 + quad1_02*l0*l2 + quad1_11*l1*l1 + quad1_12*l1*l2 + quad1_22*l2*l2;
    real quad2 = quad2_00*l0*l0 + quad2_01*l0*l1 + quad2_02*l0*l2 + quad2_11*l1*l1 + quad2_12*l1*l2 + quad2_22*l2*l2;
    real b0 = d01sq - ruc01sq - quad0;
    real b1 = d02sq - ruc02sq - quad1;
    real b2 = d03sq - ruc03sq - quad2;
    real l0_new = a00inv*b0 + a01inv*b1 + a02inv*b2;
    real l1_new = a10inv*b0 + a11inv*b1 + a12inv*b2;
    real l2_new = a20inv*b0 + a21inv*b1 + a22inv*b2;
    int almost_done = fabs(l0_new - l0) < tol && fabs(l1_new - l1) < tol && fabs(l2_new - l2) < tol;
    l0 = l0_new;
    l1 = l1_new;
    l2 = l2_new;
    if (done) break;
    done = almost_done;
  }
  l0 = l0 / dtfsq;
  l1 = l1 / dtfsq;
  l2 = l2 / dtfsq;
  cell->f[ig0] += r01 * l0;
  cell->f[ig0] += r02 * l1;
  cell->f[ig0] += r03 * l2;
  cell->f[ig1] -= r01 * l0;
  cell->f[ig2] -= r02 * l1;
  cell->f[ig3] -= r03 * l2;
}
void shake3_nodt(celldata_t *cell, int i, real tol, int maxiter){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  int ig2 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[1]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  real rmass2 = cell->rmass[ig2];
  real d01sq = cell->shake[i].rsq[0];
  real d02sq = cell->shake[i].rsq[1];
  vec<real> r01, ruc01, r02, ruc02;
  r01 = cell->x[ig0] - cell->x[ig1];
  r02 = cell->x[ig0] - cell->x[ig2];
  ruc01 = cell->x[ig0] - cell->x[ig1];
  ruc02 = cell->x[ig0] - cell->x[ig2];
  real ruc01sq = ruc01.norm();
  real ruc02sq = ruc02.norm();
  real r0101 = r01.dot(r01);
  real r0102 = r01.dot(r02);
  real r0202 = r02.dot(r02);
  real a00 = 2 * r01.dot(ruc01)*(+rmass0+rmass1);
  real a01 = 2 * r01.dot(ruc02)*(+rmass0);
  real a10 = 2 * r02.dot(ruc01)*(+rmass0);
  real a11 = 2 * r02.dot(ruc02)*(+rmass0+rmass2);
  real invmatdet = 1/(a00*a11 - a01*a10);
  real a00inv = (a11) * invmatdet;
  real a01inv = (-a01) * invmatdet;
  real a10inv = (-a10) * invmatdet;
  real a11inv = (a00) * invmatdet;
  real quad0_00 = r0101*(rmass0 + rmass1)*(rmass0 + rmass1);
  real quad0_01 = 2*r0102*rmass0*(rmass0 + rmass1);
  real quad0_11 = r0202*(rmass0)*(rmass0);
  real quad1_00 = r0101*(rmass0)*(rmass0);
  real quad1_01 = 2*r0102*rmass0*(rmass0 + rmass2);
  real quad1_11 = r0202*(rmass0 + rmass2)*(rmass0 + rmass2);
  real l0 = 0;
  real l1 = 0;
  int done = 0;
  for (int i = 0; i < maxiter; i ++) {
    real quad0 = quad0_00*l0*l0 + quad0_01*l0*l1 + quad0_11*l1*l1;
    real quad1 = quad1_00*l0*l0 + quad1_01*l0*l1 + quad1_11*l1*l1;
    real b0 = d01sq - ruc01sq - quad0;
    real b1 = d02sq - ruc02sq - quad1;
    real l0_new = a00inv*b0 + a01inv*b1;
    real l1_new = a10inv*b0 + a11inv*b1;
    int almost_done = fabs(l0_new - l0) < tol && fabs(l1_new - l1) < tol;
    l0 = l0_new;
    l1 = l1_new;
    if (done) break;
    done = almost_done;
  }
  cell->shake_xuc[ig0] += r01 * (l0 * rmass0);
  cell->shake_xuc[ig0] += r02 * (l1 * rmass0);
  cell->shake_xuc[ig1] -= r01 * (l0 * rmass1);
  cell->shake_xuc[ig2] -= r02 * (l1 * rmass2);
}
void shake3angle_nodt(celldata_t *cell, int i, real tol, int maxiter){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  int ig2 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[1]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  real rmass2 = cell->rmass[ig2];
  real d01sq = cell->shake[i].rsq[0];
  real d02sq = cell->shake[i].rsq[1];
  real d12sq = cell->shake[i].rsq[2];
  vec<real> r01, ruc01, r02, ruc02, r12, ruc12;
  r01 = cell->x[ig0] - cell->x[ig1];
  r02 = cell->x[ig0] - cell->x[ig2];
  r12 = cell->x[ig1] - cell->x[ig2];
  ruc01 = cell->x[ig0] - cell->x[ig1];
  ruc02 = cell->x[ig0] - cell->x[ig2];
  ruc12 = cell->x[ig1] - cell->x[ig2];
  real ruc01sq = ruc01.norm();
  real ruc02sq = ruc02.norm();
  real ruc12sq = ruc12.norm();
  real r0101 = r01.dot(r01);
  real r0102 = r01.dot(r02);
  real r0112 = r01.dot(r12);
  real r0202 = r02.dot(r02);
  real r0212 = r02.dot(r12);
  real r1212 = r12.dot(r12);
  real a00 = 2 * r01.dot(ruc01)*(+rmass0+rmass1);
  real a01 = 2 * r01.dot(ruc02)*(+rmass0);
  real a02 = 2 * r01.dot(ruc12)*(-rmass1);
  real a10 = 2 * r02.dot(ruc01)*(+rmass0);
  real a11 = 2 * r02.dot(ruc02)*(+rmass0+rmass2);
  real a12 = 2 * r02.dot(ruc12)*(+rmass2);
  real a20 = 2 * r12.dot(ruc01)*(-rmass1);
  real a21 = 2 * r12.dot(ruc02)*(+rmass2);
  real a22 = 2 * r12.dot(ruc12)*(+rmass1+rmass2);
  real invmatdet = 1/(a00*a11*a22 - a00*a12*a21 - a01*a10*a22 + a01*a12*a20 + a02*a10*a21 - a02*a11*a20);
  real a00inv = (a11*a22 - a12*a21) * invmatdet;
  real a01inv = (-a01*a22 + a02*a21) * invmatdet;
  real a02inv = (a01*a12 - a02*a11) * invmatdet;
  real a10inv = (-a10*a22 + a12*a20) * invmatdet;
  real a11inv = (a00*a22 - a02*a20) * invmatdet;
  real a12inv = (-a00*a12 + a02*a10) * invmatdet;
  real a20inv = (a10*a21 - a11*a20) * invmatdet;
  real a21inv = (-a00*a21 + a01*a20) * invmatdet;
  real a22inv = (a00*a11 - a01*a10) * invmatdet;
  real quad0_00 = r0101*(rmass0 + rmass1)*(rmass0 + rmass1);
  real quad0_01 = 2*r0102*rmass0*(rmass0 + rmass1);
  real quad0_02 = -2*r0112*rmass1*(rmass0 + rmass1);
  real quad0_11 = r0202*(rmass0)*(rmass0);
  real quad0_12 = -2*r0212*rmass0*rmass1;
  real quad0_22 = r1212*(rmass1)*(rmass1);
  real quad1_00 = r0101*(rmass0)*(rmass0);
  real quad1_01 = 2*r0102*rmass0*(rmass0 + rmass2);
  real quad1_02 = 2*r0112*rmass0*rmass2;
  real quad1_11 = r0202*(rmass0 + rmass2)*(rmass0 + rmass2);
  real quad1_12 = 2*r0212*rmass2*(rmass0 + rmass2);
  real quad1_22 = r1212*(rmass2)*(rmass2);
  real quad2_00 = r0101*(rmass1)*(rmass1);
  real quad2_01 = -2*r0102*rmass1*rmass2;
  real quad2_02 = -2*r0112*rmass1*(rmass1 + rmass2);
  real quad2_11 = r0202*(rmass2)*(rmass2);
  real quad2_12 = 2*r0212*rmass2*(rmass1 + rmass2);
  real quad2_22 = r1212*(rmass1 + rmass2)*(rmass1 + rmass2);
  real l0 = 0;
  real l1 = 0;
  real l2 = 0;
  int done = 0;
  for (int i = 0; i < maxiter; i ++) {
    real quad0 = quad0_00*l0*l0 + quad0_01*l0*l1 + quad0_02*l0*l2 + quad0_11*l1*l1 + quad0_12*l1*l2 + quad0_22*l2*l2;
    real quad1 = quad1_00*l0*l0 + quad1_01*l0*l1 + quad1_02*l0*l2 + quad1_11*l1*l1 + quad1_12*l1*l2 + quad1_22*l2*l2;
    real quad2 = quad2_00*l0*l0 + quad2_01*l0*l1 + quad2_02*l0*l2 + quad2_11*l1*l1 + quad2_12*l1*l2 + quad2_22*l2*l2;
    real b0 = d01sq - ruc01sq - quad0;
    real b1 = d02sq - ruc02sq - quad1;
    real b2 = d12sq - ruc12sq - quad2;
    real l0_new = a00inv*b0 + a01inv*b1 + a02inv*b2;
    real l1_new = a10inv*b0 + a11inv*b1 + a12inv*b2;
    real l2_new = a20inv*b0 + a21inv*b1 + a22inv*b2;
    int almost_done = fabs(l0_new - l0) < tol && fabs(l1_new - l1) < tol && fabs(l2_new - l2) < tol;
    l0 = l0_new;
    l1 = l1_new;
    l2 = l2_new;
    if (done) break;
    done = almost_done;
  }
  cell->shake_xuc[ig0] += r01 * (l0 * rmass0);
  cell->shake_xuc[ig0] += r02 * (l1 * rmass0);
  cell->shake_xuc[ig1] -= r01 * (l0 * rmass1);
  cell->shake_xuc[ig1] += r12 * (l2 * rmass1);
  cell->shake_xuc[ig2] -= r02 * (l1 * rmass2);
  cell->shake_xuc[ig2] -= r12 * (l2 * rmass2);
}
void shake4_nodt(celldata_t *cell, int i, real tol, int maxiter){
  int ig0 = i;
  int ig1 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[0]];
  int ig2 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[1]];
  int ig3 = cell->bonded_id[cell->first_bonded[i] + cell->shake[i].idx[2]];
  real rmass0 = cell->rmass[ig0];
  real rmass1 = cell->rmass[ig1];
  real rmass2 = cell->rmass[ig2];
  real rmass3 = cell->rmass[ig3];
  real d01sq = cell->shake[i].rsq[0];
  real d02sq = cell->shake[i].rsq[1];
  real d03sq = cell->shake[i].rsq[2];
  vec<real> r01, ruc01, r02, ruc02, r03, ruc03;
  r01 = cell->x[ig0] - cell->x[ig1];
  r02 = cell->x[ig0] - cell->x[ig2];
  r03 = cell->x[ig0] - cell->x[ig3];
  ruc01 = cell->x[ig0] - cell->x[ig1];
  ruc02 = cell->x[ig0] - cell->x[ig2];
  ruc03 = cell->x[ig0] - cell->x[ig3];
  real ruc01sq = ruc01.norm();
  real ruc02sq = ruc02.norm();
  real ruc03sq = ruc03.norm();
  real r0101 = r01.dot(r01);
  real r0102 = r01.dot(r02);
  real r0103 = r01.dot(r03);
  real r0202 = r02.dot(r02);
  real r0203 = r02.dot(r03);
  real r0303 = r03.dot(r03);
  real a00 = 2 * r01.dot(ruc01)*(+rmass0+rmass1);
  real a01 = 2 * r01.dot(ruc02)*(+rmass0);
  real a02 = 2 * r01.dot(ruc03)*(+rmass0);
  real a10 = 2 * r02.dot(ruc01)*(+rmass0);
  real a11 = 2 * r02.dot(ruc02)*(+rmass0+rmass2);
  real a12 = 2 * r02.dot(ruc03)*(+rmass0);
  real a20 = 2 * r03.dot(ruc01)*(+rmass0);
  real a21 = 2 * r03.dot(ruc02)*(+rmass0);
  real a22 = 2 * r03.dot(ruc03)*(+rmass0+rmass3);
  real invmatdet = 1/(a00*a11*a22 - a00*a12*a21 - a01*a10*a22 + a01*a12*a20 + a02*a10*a21 - a02*a11*a20);
  real a00inv = (a11*a22 - a12*a21) * invmatdet;
  real a01inv = (-a01*a22 + a02*a21) * invmatdet;
  real a02inv = (a01*a12 - a02*a11) * invmatdet;
  real a10inv = (-a10*a22 + a12*a20) * invmatdet;
  real a11inv = (a00*a22 - a02*a20) * invmatdet;
  real a12inv = (-a00*a12 + a02*a10) * invmatdet;
  real a20inv = (a10*a21 - a11*a20) * invmatdet;
  real a21inv = (-a00*a21 + a01*a20) * invmatdet;
  real a22inv = (a00*a11 - a01*a10) * invmatdet;
  real quad0_00 = r0101*(rmass0 + rmass1)*(rmass0 + rmass1);
  real quad0_01 = 2*r0102*rmass0*(rmass0 + rmass1);
  real quad0_02 = 2*r0103*rmass0*(rmass0 + rmass1);
  real quad0_11 = r0202*(rmass0)*(rmass0);
  real quad0_12 = 2*r0203*(rmass0)*(rmass0);
  real quad0_22 = r0303*(rmass0)*(rmass0);
  real quad1_00 = r0101*(rmass0)*(rmass0);
  real quad1_01 = 2*r0102*rmass0*(rmass0 + rmass2);
  real quad1_02 = 2*r0103*(rmass0)*(rmass0);
  real quad1_11 = r0202*(rmass0 + rmass2)*(rmass0 + rmass2);
  real quad1_12 = 2*r0203*rmass0*(rmass0 + rmass2);
  real quad1_22 = r0303*(rmass0)*(rmass0);
  real quad2_00 = r0101*(rmass0)*(rmass0);
  real quad2_01 = 2*r0102*(rmass0)*(rmass0);
  real quad2_02 = 2*r0103*rmass0*(rmass0 + rmass3);
  real quad2_11 = r0202*(rmass0)*(rmass0);
  real quad2_12 = 2*r0203*rmass0*(rmass0 + rmass3);
  real quad2_22 = r0303*(rmass0 + rmass3)*(rmass0 + rmass3);
  real l0 = 0;
  real l1 = 0;
  real l2 = 0;
  int done = 0;
  for (int i = 0; i < maxiter; i ++) {
    real quad0 = quad0_00*l0*l0 + quad0_01*l0*l1 + quad0_02*l0*l2 + quad0_11*l1*l1 + quad0_12*l1*l2 + quad0_22*l2*l2;
    real quad1 = quad1_00*l0*l0 + quad1_01*l0*l1 + quad1_02*l0*l2 + quad1_11*l1*l1 + quad1_12*l1*l2 + quad1_22*l2*l2;
    real quad2 = quad2_00*l0*l0 + quad2_01*l0*l1 + quad2_02*l0*l2 + quad2_11*l1*l1 + quad2_12*l1*l2 + quad2_22*l2*l2;
    real b0 = d01sq - ruc01sq - quad0;
    real b1 = d02sq - ruc02sq - quad1;
    real b2 = d03sq - ruc03sq - quad2;
    real l0_new = a00inv*b0 + a01inv*b1 + a02inv*b2;
    real l1_new = a10inv*b0 + a11inv*b1 + a12inv*b2;
    real l2_new = a20inv*b0 + a21inv*b1 + a22inv*b2;
    int almost_done = fabs(l0_new - l0) < tol && fabs(l1_new - l1) < tol && fabs(l2_new - l2) < tol;
    l0 = l0_new;
    l1 = l1_new;
    l2 = l2_new;
    if (done) break;
    done = almost_done;
  }
  cell->shake_xuc[ig0] += r01 * (l0 * rmass0);
  cell->shake_xuc[ig0] += r02 * (l1 * rmass0);
  cell->shake_xuc[ig0] += r03 * (l2 * rmass0);
  cell->shake_xuc[ig1] -= r01 * (l0 * rmass1);
  cell->shake_xuc[ig2] -= r02 * (l1 * rmass2);
  cell->shake_xuc[ig3] -= r03 * (l2 * rmass3);
}
